Skip to content

Examples and Plotting

Plots, simple¤

Plots.jl and Makie.jl are fully supported by Rasters.jl, with recipes for plotting Raster and RasterStack provided. plot will plot a heatmap with axes matching dimension values. If mappedcrs is used, converted values will be shown on axes instead of the underlying crs values. contourf will similarly plot a filled contour plot.

Pixel resolution is limited to allow loading very large files quickly. max_res specifies the maximum pixel resolution to show on the longest axis of the array. It can be set manually to change the resolution (e.g. for large or high-quality plots):

using Rasters, RasterDataSources, ArchGDAL, Plots
A = Raster(WorldClim{BioClim}, 5)
plot(A; max_res=3000)

For Makie, plot functions in a similar way. plot will only accept two-dimensional rasters. You can invoke contour, contourf, heatmap, surface or any Makie plotting function which supports surface-like data on a 2D raster.

To obtain tiled plots for 3D rasters and RasterStacks, use the function Rasters.rplot([gridposition], raster; kw_args...). This is an unexported function, since we're not sure how the API will change going forward.

Makie, simple¤

using CairoMakie
CairoMakie.activate!(px_per_unit = 2)
using Rasters, CairoMakie, RasterDataSources, ArchGDAL
A = Raster(WorldClim{BioClim}, 5)
Makie.plot(A)

Loading data¤

Our first example simply loads a file from disk and plots it.

This netcdf file only has one layer, if it has more we could use RasterStack instead.

using Rasters, NCDatasets, Plots
url = "https://www.unidata.ucar.edu/software/netcdf/examples/tos_O1_2001-2002.nc";
filename = download(url, "tos_O1_2001-2002.nc");
A = Raster(filename)

180×170×24 Raster{Union{Missing, Float32},3} tos with dimensions: 
  X Mapped{Float64} Float64[1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
  Y Mapped{Float64} Float64[-79.5, -78.5, …, 88.5, 89.5] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
  Ti Sampled{DateTime360Day} DateTime360Day[DateTime360Day(2001-01-16T00:00:00), …, DateTime360Day(2002-12-16T00:00:00)] ForwardOrdered Explicit Intervals
extent: Extent(X = (-0.0, 360.0), Y = (-80.0, 90.0), Ti = (DateTime360Day(2001-01-01T00:00:00), DateTime360Day(2003-01-01T00:00:00)))missingval: missingcrs: EPSG:4326
mappedcrs: EPSG:4326
parent:
[:, :, 1]
        -79.5       -78.5       …   86.5     87.5     88.5     89.5
   1.0     missing     missing     271.43   271.437  271.445  271.459
   3.0     missing     missing     271.431  271.438  271.445  271.459
   5.0     missing     missing     271.431  271.438  271.445  271.459
   7.0     missing     missing     271.431  271.439  271.446  271.459
   ⋮                            ⋱                      ⋮      
 351.0     missing     missing     271.43   271.435  271.445  271.459
 353.0     missing     missing     271.43   271.436  271.445  271.459
 355.0     missing     missing     271.43   271.436  271.445  271.459
 357.0     missing     missing     271.43   271.437  271.445  271.459
 359.0     missing     missing  …  271.43   271.437  271.445  271.459
[and 23 more slices...]

Objects with Dimensions other than X and Y will produce multi-pane plots. Here we plot every third month in the first year in one plot:

A[Ti=1:3:12] |> plot

Now plot the ocean temperatures around the Americas in the first month of 2001. Notice we are using lat/lon coordinates and date/time instead of regular indices. The time dimension uses DateTime360Day, so we need to load CFTime.jl to index it with Near.

using CFTime
A[Ti(Near(DateTime360Day(2001, 01, 17))), Y(-60.0 .. 90.0), X(45.0 .. 190.0)] |> plot

Now get the mean over the timespan, then save it to disk, and plot it as a filled contour.

Other plot functions and sliced objects that have only one X/Y/Z dimension fall back to generic DimensionalData.jl plotting, which will still correctly label plot axes.

using Statistics
# Take the mean
mean_tos = mean(A; dims=Ti)

180×170×1 Raster{Union{Missing, Float32},3} tos with dimensions: 
  X Mapped{Float64} Float64[1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
  Y Mapped{Float64} Float64[-79.5, -78.5, …, 88.5, 89.5] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
  Ti Sampled{DateTime360Day} DateTime360Day(2002-01-16T00:00:00):Millisecond(2592000000):DateTime360Day(2002-01-16T00:00:00) ForwardOrdered Explicit Intervals
extent: Extent(X = (-0.0, 360.0), Y = (-80.0, 90.0), Ti = (DateTime360Day(2001-01-01T00:00:00), DateTime360Day(2003-01-01T00:00:00)))missingval: missingcrs: EPSG:4326
mappedcrs: EPSG:4326
parent:
[:, :, 1]
        -79.5       -78.5       …   86.5     87.5     88.5     89.5
   1.0     missing     missing     271.427  271.434  271.443  271.454
   3.0     missing     missing     271.427  271.434  271.443  271.454
   5.0     missing     missing     271.427  271.434  271.443  271.454
   7.0     missing     missing     271.427  271.435  271.444  271.454
   ⋮                            ⋱                      ⋮      
 351.0     missing     missing     271.427  271.433  271.443  271.454
 353.0     missing     missing     271.427  271.433  271.443  271.454
 355.0     missing     missing     271.427  271.433  271.443  271.454
 357.0     missing     missing     271.427  271.433  271.443  271.454
 359.0     missing     missing  …  271.427  271.433  271.443  271.454

Plot a contour plot¤

using Plots
Plots.contourf(mean_tos; dpi=300, size=(800, 400))

write to disk¤

Write the mean values to disk

write("mean_tos.nc", mean_tos)

"mean_tos.nc"

Plotting recipes in DimensionalData.jl are the fallback for Rasters.jl when the object doesn't have 2 X/Y/Z dimensions, or a non-spatial plot command is used. So (as a random example) we could plot a transect of ocean surface temperature at 20 degree latitude :

A[Y(Near(20.0)), Ti(1)] |> plot

Polygon masking, mosaic and plot¤

In this example we will mask the Scandinavian countries with border polygons, then mosaic together to make a single plot.

First, get the country boundary shape files using GADM.jl.

using Rasters, RasterDataSources, ArchGDAL, Shapefile, Plots, Dates, Downloads, NCDatasets

WARNING: using Downloads.download in module Main conflicts with an existing identifier.

Download the shapefile

shapefile_url = "https://github.com/nvkelso/natural-earth-vector/raw/master/10m_cultural/ne_10m_admin_0_countries.shp"
shapefile_name = "boundary_lines.shp"
Downloads.download(shapefile_url, shapefile_name)

"boundary_lines.shp"

Load using Shapefile.jl

shapes = Shapefile.Handle(shapefile_name)
denmark_border = shapes.shapes[71]
norway_border = shapes.shapes[53]
sweden_border = shapes.shapes[54]

Polygon(4665 Points)

Then load raster data. We load some worldclim layers using RasterDataSources via Rasters.jl:

climate = RasterStack(WorldClim{Climate}, (:tmin, :tmax, :prec, :wind); month=July)

RasterStack with dimensions: 
  X Projected{Float64} LinRange{Float64}(-180.0, 179.833, 2160) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(89.8333, -90.0, 1080) ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
  :tmin Float32 dims: X, Y (2160×1080)
  :tmax Float32 dims: X, Y (2160×1080)
  :prec Int16 dims: X, Y (2160×1080)
  :wind Float32 dims: X, Y (2160×1080)

mask Denmark, Norway and Sweden from the global dataset using their border polygon, then trim the missing values. We pad trim with a 10 pixel margin.

mask_trim(climate, poly) = trim(mask(climate; with=poly); pad=10)

denmark = mask_trim(climate, denmark_border)
norway = mask_trim(climate, norway_border)
sweden = mask_trim(climate, sweden_border)

RasterStack with dimensions: 
  X Projected{Float64} LinRange{Float64}(9.5, 25.6667, 98) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(70.5, 53.6667, 102) ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
  :tmin Float32 dims: X, Y (98×102)
  :tmax Float32 dims: X, Y (98×102)
  :prec Int16 dims: X, Y (98×102)
  :wind Float32 dims: X, Y (98×102)

Then load raster data. We load some worldclim layers using RasterDataSources via Rasters.jl:

climate = RasterStack(WorldClim{Climate}, (:tmin, :tmax, :prec, :wind); month=July)

RasterStack with dimensions: 
  X Projected{Float64} LinRange{Float64}(-180.0, 179.833, 2160) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(89.8333, -90.0, 1080) ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
  :tmin Float32 dims: X, Y (2160×1080)
  :tmax Float32 dims: X, Y (2160×1080)
  :prec Int16 dims: X, Y (2160×1080)
  :wind Float32 dims: X, Y (2160×1080)

mask Denmark, Norway and Sweden from the global dataset using their border polygon, then trim the missing values. We pad trim with a 10 pixel margin.

mask_trim(climate, poly) = trim(mask(climate; with=poly); pad=10)

denmark = mask_trim(climate, denmark_border)
norway = mask_trim(climate, norway_border)
sweden = mask_trim(climate, sweden_border)

RasterStack with dimensions: 
  X Projected{Float64} LinRange{Float64}(9.5, 25.6667, 98) ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(70.5, 53.6667, 102) ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
  :tmin Float32 dims: X, Y (98×102)
  :tmax Float32 dims: X, Y (98×102)
  :prec Int16 dims: X, Y (98×102)
  :wind Float32 dims: X, Y (98×102)

Plotting with Plots¤

First define a function to add borders to all subplots.

function borders!(p, poly)
    for i in 1:length(p)
        Plots.plot!(p, poly; subplot=i, fillalpha=0, linewidth=0.6)
    end
    return p
end

borders! (generic function with 1 method)

Now we can plot the individual countries.

dp = plot(denmark)
borders!(dp, denmark_border)

and sweden

sp = plot(sweden)
borders!(sp, sweden_border)

and norway

np = plot(norway)
borders!(np, norway_border)

The Norway shape includes a lot of islands. Lets crop them out using .. intervals:

norway_region = climate[X(0..40), Y(55..73)]
plot(norway_region)

And mask it with the border again:

norway = mask_trim(norway_region, norway_border)
np = plot(norway)
borders!(np, norway_border)

Now we can combine the countries into a single raster using mosaic. first will take the first value if/when there is an overlap.

scandinavia = mosaic(first, denmark, norway, sweden)

RasterStack with dimensions: 
  X Projected{Float64} 3.1666666666666443:0.16666666666666666:32.49999999999998 ForwardOrdered Regular Intervals crs: WellKnownText,
  Y Projected{Float64} 72.66666666666666:-0.16666666666666666:52.99999999999999 ReverseOrdered Regular Intervals crs: WellKnownText
and 4 layers:
  :tmin Float32 dims: X, Y (177×119)
  :tmax Float32 dims: X, Y (177×119)
  :prec Int16 dims: X, Y (177×119)
  :wind Float32 dims: X, Y (177×119)

And plot scandinavia, with all borders included:

p = plot(scandinavia)
borders!(p, denmark_border)
borders!(p, norway_border)
borders!(p, sweden_border)
p

And save to netcdf - a single multi-layered file, and tif, which will write a file for each stack layer.

write("scandinavia.nc", scandinavia)
write("scandinavia.tif", scandinavia)

(tmin = "scandinavia_tmin.tif", tmax = "scandinavia_tmax.tif", prec = "scandinavia_prec.tif", wind = "scandinavia_wind.tif")

Rasters.jl provides a range of other methods that are being added to over time. Where applicable these methods read and write lazily to and from disk-based arrays of common raster file types. These methods also work for entire RasterStacks and RasterSeries using the same syntax.

This page was generated using Literate.jl.